% This file contains code for empirical / graphical component of the
% Handel, Kolstad, Spinniwijn paper on Adverse Selection and Information Frictions in Insurance Markets. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Section 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define survey vars
SvyDummies = cat(5, BenefitsKnowledgeAnyIncorrect,BenefitsKnowledgeAnyNotSure,TimeCost,TimeCost.*TimeCostAccept, TimeCost.*TimeCostDontLike,ProviderNetworkHSPMore, ProviderNetworkPPOMore, ProviderNetworkNotSure,TMEGuessOverestimate, TMEGuessUnderestimate, TMEGuessNotSure,TaxBenefitsMisunderstands,TaxBenefitsNotSure);

% Desubsidize the HSP
MSsubsidy = zeros(size(UtilityHSP));
index1 = find(familysizes==1);
index2 = find(familysizes==2);
index3 = find(min(familysizes,3)==3);
MSsubsidy(index1) = 1500;
MSsubsidy(index2) = 3000;
MSsubsidy(index3) = 3750;

HSP = UtilityHSP - MSsubsidy;

% People who chose PPO
finder = find(choices(:,2)==0);
varlist = {'RA','E','HSP','UtilityPPO','ages','incomes','plans','genders','familysizes','famdraws','famdrawsOOP'};
for i = 1:length(varlist)
    eval([varlist{i},' = ',varlist{i},'(finder,2,:,:);'])
end
SvyDummies = SvyDummies(finder,2,:,:,:);

% Use parameters from Handel and Kolstad (2015) "Health Insurance for "Humans"" - Full Model with Frictions and Inertia
Params = FullModelwHSAAlpha;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%  Figures  %%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Use tier 3 for analysis
finder = find(min(familysizes(:,1,1,1),3)==3);
alphalist = 0:0.1:1;

% WTP to buy the PPO above the HSP
WTPs = cell(numel(alphalist),1);
for alphaIND = 1:numel(alphalist)
	alpha = alphalist(alphaIND);
	% Calculate WTP under adjusted frictions
	NewParams = Params;
	NewParams(7:numel(NewParams)) = (1-alpha) * NewParams(7:numel(NewParams));
	ComputedWTP = 0 - CertaintyEquivalence(NewParams,RA,E,HSP,ages,genders,incomes,familysizes,plans,SvyDummies,[],'true');
	WTPs{alphaIND} = ComputedWTP(finder,:,:);
end
OOP = famdrawsOOP(finder,:,:,:);
allowed = famdraws(finder,:,:,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Figures: Mkt Equilibrium with different friction levels %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for alpha = [0,0.5,1]
	alphaIND = find(alphalist==alpha);
	if alpha==1
		alphaSTR = '10';
	elseif alpha==0.5
		alphaSTR = '05';
	else
		alphaSTR = '00';
	end
	
	WTP = WTPs{alphaIND};
	FriclessWTP = WTPs{find(alphalist==1)};
	
	% Demand = WTP for PPO
	[Demand,ind] = sort(WTP(:),'descend');

	% Marginal cost = expected cost of the next buyer, where next is ordered by WTP.
	meanOOP = mean(OOP,3);
	MarginalCost = meanOOP(ind);
    
	% Avg cost = sum of marginal costs, over the # of buyers
	AverageCost = cumsum(MarginalCost) ./ (1:length(MarginalCost))';

    % Value is demand without frictions. Plotting conditional on the employee's position on the with-frictions demand curve.
    Value = FriclessWTP(ind);

    % Figure - Etc.
	xaxis = (1:length(MarginalCost))' / length(MarginalCost);
	
            %%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%% Spline plot %%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%
	
	% Define the number of splines, then divide the x-axis accordingly
	numsplines = 10;
	
    % 11 groups, where the first and last are half the size of the other 9 groups, and form the intercept and the end point
	splinexaxis = floor((xaxis*numsplines) + .5 - 10^-6) + 1;
	
	% Mean spending/demand within each group
	DemandSpline = accumarray(splinexaxis,Demand,[],@mean);
	AverageCostSpline = accumarray(splinexaxis,AverageCost,[],@mean);
	MarginalCostSpline = accumarray(splinexaxis,MarginalCost,[],@mean);
	
    % Adjust so AC and MC start at the same intercept, same w D and MR
	MarginalCostSpline = MarginalCostSpline + (AverageCostSpline(1) - MarginalCostSpline(1));
	
    % Compute avg. value premium above MC for each group, and add it to MC
	RiskPremium = accumarray(splinexaxis,Value - MarginalCost,[],@mean);
	ValueSpline = MarginalCostSpline + RiskPremium;

    % Figure - Etc.
	splinexaxis = ((accumarray(splinexaxis,splinexaxis,[],@max)) - 1) / numsplines;
	if alpha~=1
		plot(splinexaxis,DemandSpline,'k-',splinexaxis,AverageCostSpline,'k--',splinexaxis,MarginalCostSpline,'k-.',splinexaxis,ValueSpline,'k:')
		legend('Demand','Average Cost','Marginal Cost','Value')
	else
		plot(splinexaxis,DemandSpline,'k-',splinexaxis,AverageCostSpline,'k--',splinexaxis,MarginalCostSpline,'k-.')
		legend('Demand','Average Cost','Marginal Cost')
    end
	xlabel('Quantity')
	ylabel('Value')
	axis([0 1 2000 12000])
	outfile = ['MktEquilibriumAlpha',alphaSTR,'FixedAxes.png'];
	saveas(figure(1),outfile)
		
	eval(['Alpha',num2str(alpha*10),'AverageCostSpline = AverageCostSpline;'])
end
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Figure: Avg Cost curves with different friction levels %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot(splinexaxis,Alpha0AverageCostSpline,'k-',splinexaxis,Alpha5AverageCostSpline,'k--',splinexaxis,Alpha10AverageCostSpline,'k-.')
legend('\alpha = 0','\alpha = 0.5','\alpha = 1')
xlabel('Quantity')
ylabel('Cost')

outfile = ['AvgCostCurves.png'];
saveas(figure(1),outfile)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Figures: Mkt equilibrium with risk adjustment %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for alpha = [0,1]
	if alpha==0
		alphaSTR = '00';
	else
		alphaSTR = '10';
	end
	
	WTP = WTPs{find(alphalist==alpha)};
	
	% Get Demand = WTP for PPO, sort descending
	[Demand,ind] = sort(WTP(:),'descend');
	AverageCostRiskAdjusted = cell(3,1);
	for betaIND = 1:3
		beta = (betaIND - 1)/2;
        
		% Marginal cost = expected cost of the next buyer, where next is ordered by WTP.
		meanOOP = (1-beta)*mean(OOP,3) + beta*mean(OOP(:));
		MarginalCost = meanOOP(ind);
		
        % Avg cost = sum of marginal costs, over the # of buyers
		AverageCostRiskAdjusted{betaIND} = cumsum(MarginalCost) ./ (1:length(MarginalCost))';
    end
		splinexaxis = floor((xaxis*numsplines) + .5 - 10^-6) + 1;
	
	%Get the mean spending/demand within each group
	DemandSpline = accumarray(splinexaxis,Demand,[],@mean);
	AverageCostRiskAdjustedSpline = cell(3,1);
	for betaIND = 1:3
		AverageCostRiskAdjustedSpline{betaIND} = accumarray(splinexaxis,AverageCostRiskAdjusted{betaIND},[],@mean);
    end
	splinexaxis = ((accumarray(splinexaxis,splinexaxis,[],@max)) - 1) / numsplines;
	plot(splinexaxis,DemandSpline,'k-',splinexaxis,AverageCostRiskAdjustedSpline{1},'k--',splinexaxis,AverageCostRiskAdjustedSpline{2},'k-.',splinexaxis,AverageCostRiskAdjustedSpline{3},'k:')
	legend('Demand','Average Cost, \beta = 0','Average Cost, \beta = 0.5','Average Cost, \beta = 1')

	outfile = ['RiskAdjustmentMktEquilibriumAlpha',alphaSTR,'.png'];
	saveas(figure(1),outfile)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 3D plots of outcomes under different policies: PPO share premium, welfare %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Shares = NaN(11,11);
Premiums = NaN(11,11);
Welfares = NaN(11,11);

for alphaIND = 1:11
for betaIND = 1:11
	alpha = (alphaIND - 1) / 10;	
	beta = (betaIND - 1) / 10;
	WTP = WTPs{alphaIND};
	
    %Get Demand = WTP for PPO, sort descending
	[Demand,ind] = sort(WTP(:),'descend');

	% Marginal cost = expected cost of the next buyer, where next is ordered by WTP.
	meanOOP = (1-beta)*mean(OOP,3) + beta*mean(OOP(:));
	MarginalCost = meanOOP(ind);

    % Avg cost = sum of marginal costs, over the # of buyers
	AverageCost = cumsum(MarginalCost) ./ (1:length(MarginalCost))';
	
	%Calculate unadjusted MC for use in welfare
	meanOOPTrue = mean(OOP,3) ;
	MarginalCostTrue = meanOOPTrue(ind);
	Value = WTPs{find(alphalist==1)};
	Value = Value(ind);
	Share = mean(Demand >= AverageCost);
	Premium = min(AverageCost(Demand >= AverageCost));
	Welfare = mean((Value - MarginalCostTrue).*(Demand >= AverageCost));
	Shares(alphaIND,betaIND) = Share;
	Premiums(alphaIND,betaIND) = Premium;
	Welfares(alphaIND,betaIND) = Welfare;
end
end

% Welfare impact relative to (0,0)
WelfareImpacts = Welfares - Welfares(1,1);

% 3D Plots
surf([0:.1:1],[0:.1:1],Shares')
colormap(gray)
xlabel('Alpha')
ylabel('Beta')
zlabel('Share')
saveas(figure(1),'3DGridShares.png')
surf([0:.1:1],[0:.1:1],Premiums')
colormap(gray)
xlabel('Alpha')
ylabel('Beta')
zlabel('Premium Difference')
saveas(figure(1),'3DGridPremiums.png')
surf([0:.1:1],[0:.1:1],WelfareImpacts')
colormap(gray)
xlabel('Alpha')
ylabel('Beta')
zlabel('Welfare Impact')
saveas(figure(1),'3DGridWelfareImpacts.png')

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%% Add. stat in different policy equilibria %%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% How many people are making mistakes?
UnderInsured = NaN(11,11);
OverInsured = NaN(11,11);
MeanCS = NaN(11,11);
StDevCS = NaN(11,11);
MeanCSBuyers = NaN(11,11);
StDevCSBuyers = NaN(11,11);
MeanCSNonbuyers = NaN(11,11);
StDevCSNonbuyers = NaN(11,11);
MeanWelfareBuyers = NaN(11,11);
StDevWelfareBuyers = NaN(11,11);
MeanWelfareNonbuyers = NaN(11,11);
StDevWelfareNonbuyers = NaN(11,11);

for alphaIND = 1:11
for betaIND = 1:11
	alphaSTR = num2str(alphaIND - 1);
	betaSTR = num2str(betaIND - 1);
	if alphaIND < 11
		alphaSTR = ['0',alphaSTR];
	end
	if betaIND < 11
		betaSTR = ['0',betaSTR];
	end
	
	Price = Premiums(alphaIND,betaIND);
	WTP = WTPs{alphaIND};
	WTP = WTP(:);
	Value = WTPs{find(alphalist==1)};
	Value = Value(:);
	MarginalCost = (1-beta)*mean(OOP,3) + beta*mean(OOP(:));
	MarginalCost = MarginalCost(:);
	
	UnderInsured(alphaIND,betaIND) = mean((WTP < Price).*(Value >= Price));
	OverInsured(alphaIND,betaIND) = mean((WTP >= Price).*(Value <= Price));
	
	MeanCS(alphaIND,betaIND) = mean(Value - Price);
	StDevCS(alphaIND,betaIND) = std(Value - Price);
	ksdensity(Value - Price)
	xlabel('Consumer Surplus')
	ylabel('Density')
	saveas(figure(1),['DensityCSAlpha',alphaSTR,'Beta',betaSTR,'.png'])
	
	MeanCSBuyers(alphaIND,betaIND) = mean(Value(WTP >= Price) - Price);
	StDevCSBuyers(alphaIND,betaIND) = std(Value(WTP >= Price) - Price);
	ksdensity(Value(WTP >= Price) - Price)
	xlabel('Consumer Surplus')
	ylabel('Density')
	saveas(figure(1),['DensityCSBuyersAlpha',alphaSTR,'Beta',betaSTR,'.png'])
	
	MeanCSNonbuyers(alphaIND,betaIND) = mean(Value(WTP < Price) - Price);
	StDevCSNonbuyers(alphaIND,betaIND) = std(Value(WTP < Price) - Price);
	ksdensity(Value(WTP < Price) - Price)
	xlabel('Consumer Surplus')
	ylabel('Density')
	saveas(figure(1),['DensityCSNonbuyersAlpha',alphaSTR,'Beta',betaSTR,'.png'])
	
	MeanWelfareBuyers(alphaIND,betaIND) = mean(Value(WTP >= Price) - MarginalCost(WTP >= Price));
	StDevWelfareBuyers(alphaIND,betaIND) = std(Value(WTP >= Price) - MarginalCost(WTP >= Price));
	ksdensity(Value(WTP >= Price) - MarginalCost(WTP >= Price))
	xlabel('Welfare')
	ylabel('Density')
	saveas(figure(1),['DensityWelfareBuyersAlpha',alphaSTR,'Beta',betaSTR,'.png'])
	
	MeanWelfareNonbuyers(alphaIND,betaIND) = mean(Value(WTP < Price) - MarginalCost(WTP < Price));
	StDevWelfareNonbuyers(alphaIND,betaIND) = std(Value(WTP < Price) - MarginalCost(WTP < Price));
	ksdensity(Value(WTP < Price) - MarginalCost(WTP < Price))
	xlabel('Welfare')
	ylabel('Density')
	saveas(figure(1),['DensityWelfareNonbuyersAlpha',alphaSTR,'Beta',betaSTR,'.png'])
end
end

            %%%%%%%%%%%%%%%%%%%%%%%
            %%%  Miscellaneous  %%%
            %%%%%%%%%%%%%%%%%%%%%%%
            
%%%% Distribution of PPO consumer surplus from risk protection

% WTP for full model
WTP = 0 - CertaintyEquivalence(Params,RA,E,HSP,ages,genders,incomes,familysizes,plans,SvyDummies,[],'true');
% Value for full model
Value = 0 - CertaintyEquivalence(Params(1:6),RA,E,HSP,ages,genders,incomes,familysizes,plans,[],[],'true');
% WTP/value (same thing) for naive model
ValueNaive = 0 - CertaintyEquivalence(NaiveModelwHSAAlpha,RA,E,HSP,ages,genders,incomes,familysizes,plans,[],[],'true');
% WTP/value (same thing) for extreme naive model
ValueExtremeNaive = 0 - CertaintyEquivalence(NaiveModelNoInertiawHSAAlpha,RA,E,HSP,ages,genders,incomes,familysizes,plans,[],[],'true');
% MC for insurer
MarginalCost = mean(famdrawsOOP,3);
%Expected costs for individual
ExpectedCost = mean(famdraws,3);

% Sum up frictions
Frictions = zeros(size(UtilityPPO));
SurveyParams = Params(7:length(Params));
if ~isempty(SvyDummies)
	for i = 1:size(SvyDummies,5)
		Frictions = Frictions + (SurveyParams(i)*1000*SvyDummies(:,:,:,:,i));
	end
end
Frictions = mean(Frictions,3);

WTPCell = cell(3,1);
ValueCell = cell(3,1);
ValueNaiveCell = cell(3,1);
ValueExtremeNaiveCell = cell(3,1);
MarginalCostCell = cell(3,1);
ExpectedCostCell = cell(3,1);
WelfareCell = cell(3,1);
WelfareExtremeNaiveCell = cell(3,1);
FrictionsCell = cell(3,1);
for tier = 1:3
	finder = find(min(familysizes(:,1,1,1),3)==tier);
	WTPTier = WTP(finder,:,:,:);
	ValueTier = Value(finder,:,:,:);
	ValueNaiveTier = ValueNaive(finder,:,:,:);
	ValueExtremeNaiveTier = ValueExtremeNaive(finder,:,:,:);
	MarginalCostTier = MarginalCost(finder,:,:,:);
	ExpectedCostTier = ExpectedCost(finder,:,:,:);
	FrictionsTier = Frictions(finder,:,:,:);
	
	WTPCell{tier} = WTPTier(:);
	ValueCell{tier} = ValueTier(:);
	ValueNaiveCell{tier} = ValueNaiveTier(:);
	ValueExtremeNaiveCell{tier} = ValueExtremeNaiveTier(:);
	MarginalCostCell{tier} = MarginalCostTier(:);
	ExpectedCostCell{tier} = ExpectedCostTier(:);
	FrictionsCell{tier} = FrictionsTier(:);
	WelfareCell{tier} = ValueTier(:) - MarginalCostTier(:);
	WelfareExtremeNaiveCell{tier} = ValueNaiveTier(:) - MarginalCostTier(:);
end

statfile = fopen(DistributionalStatistics.txt','w');

% Write means and std deviations into file
fprintf(statfile,'WTP for PPO, full model\r\n');
for tier = 1:3
	fprintf(statfile,'Tier %i Mean = %4.2f\r\nTier %i Std. Dev. = %4.2f\r\n',tier,mean(WTPCell{tier}),tier,std(WTPCell{tier}));
end
fprintf(statfile,'PPO value, full model\r\n');
for tier = 1:3
	fprintf(statfile,'Tier %i Mean = %4.2f\r\nTier %i Std. Dev. = %4.2f\r\n',tier,mean(ValueCell{tier}),tier,std(ValueCell{tier}));
end
fprintf(statfile,'WTP for PPO/PPO value, naive model\r\n');
for tier = 1:3
	fprintf(statfile,'Tier %i Mean = %4.2f\r\nTier %i Std. Dev. = %4.2f\r\n',tier,mean(ValueNaiveCell{tier}),tier,std(ValueNaiveCell{tier}));
end
fprintf(statfile,'WTP for PPO/PPO value, extreme naive model\r\n');
for tier = 1:3
	fprintf(statfile,'Tier %i Mean = %4.2f\r\nTier %i Std. Dev. = %4.2f\r\n',tier,mean(ValueExtremeNaiveCell{tier}),tier,std(ValueExtremeNaiveCell{tier}));
end
fprintf(statfile,'Marginal costs\r\n');
for tier = 1:3
	fprintf(statfile,'Tier %i Mean = %4.2f\r\nTier %i Std. Dev. = %4.2f\r\n',tier,mean(MarginalCostCell{tier}),tier,std(MarginalCostCell{tier}));
end
fprintf(statfile,'Expected costs\r\n');
for tier = 1:3
	fprintf(statfile,'Tier %i Mean = %4.2f\r\nTier %i Std. Dev. = %4.2f\r\n',tier,mean(ExpectedCostCell{tier}),tier,std(ExpectedCostCell{tier}));
end
fprintf(statfile,'PPO welfare, full model\r\n');
for tier = 1:3
	fprintf(statfile,'Tier %i Mean = %4.2f\r\nTier %i Std. Dev. = %4.2f\r\n',tier,mean(WelfareCell{tier}),tier,std(WelfareCell{tier}));
end
fprintf(statfile,'PPO welfare, naive model\r\n');
for tier = 1:3
	fprintf(statfile,'Tier %i Mean = %4.2f\r\nTier %i Std. Dev. = %4.2f\r\n',tier,mean(WelfareExtremeNaiveCell{tier}),tier,std(WelfareExtremeNaiveCell{tier}));
end
fprintf(statfile,'Frictions\r\n');
for tier = 1:3
	fprintf(statfile,'Tier %i Mean = %4.2f\r\nTier %i Std. Dev. = %4.2f\r\n',tier,mean(FrictionsCell{tier}),tier,std(FrictionsCell{tier}));
end
fclose(statfile);

% Density plots
for tier = 1:3
	ksdensity(WTPCell{tier})
	xlabel('Value')
	ylabel('Density')
	saveas(figure(1),['WTPDensityTier',num2str(tier),'.png'])
end
for tier = 1:3
	ksdensity(ValueCell{tier})
	xlabel('Value')
	ylabel('Density')
	saveas(figure(1),['ValueDensityTier',num2str(tier),'.png'])
end
for tier = 1:3
	ksdensity(ValueNaiveCell{tier})
	xlabel('Value')
	ylabel('Density')
	saveas(figure(1),['ValueNaiveDensityTier',num2str(tier),'.png'])
end
for tier = 1:3
	ksdensity(ValueExtremeNaiveCell{tier})
	xlabel('Value')
	ylabel('Density')
	saveas(figure(1),['ValueExtremeNaiveDensityTier',num2str(tier),'.png'])
end
for tier = 1:3
	ksdensity(MarginalCostCell{tier})
	xlabel('Cost')
	ylabel('Density')
	saveas(figure(1),['MarginalCostTier',num2str(tier),'.png'])
end
for tier = 1:3
	ksdensity(ExpectedCostCell{tier})
	xlabel('Cost')
	ylabel('Density')
	saveas(figure(1),['ExpectedCostTier',num2str(tier),'.png'])
end
for tier = 1:3
	ksdensity(WelfareCell{tier})
	xlabel('Welfare')
	ylabel('Density')
	saveas(figure(1),['WelfareDensityTier',num2str(tier),'.png'])
end
for tier = 1:3
	ksdensity(WelfareExtremeNaiveCell{tier})
	xlabel('Welfare')
	ylabel('Density')
	saveas(figure(1),['WelfareExtremeNaiveDensityTier',num2str(tier),'.png'])
end
for tier = 1:3
	ksdensity(FrictionsCell{tier})
	xlabel('Friction')
	ylabel('Density')
	saveas(figure(1),['FrictionsDensityTier',num2str(tier),'.png'])
end

% Correlations between variables
CorrelationsCell = cell(3,1);
for tier = 1:3
    Corr = corr([WTPCell{tier},ValueCell{tier},ValueNaiveCell{tier},ValueExtremeNaiveCell{tier},MarginalCostCell{tier},ExpectedCostCell{tier},WelfareCell{tier},WelfareExtremeNaiveCell{tier},FrictionsCell{tier}]);
	CorrelationsCell{tier} = Corr;
end
save(Correlations.mat','CorrelationsCell')

% "Marginal effects"
WTP = WTPs{find(alphalist==0)};
Value = WTPs{find(alphalist==1)};
offset = 10;

ksdensity(Value(WTP(:) >= prctile(WTP(:),50-offset) & WTP(:) <= prctile(WTP(:),50+offset)))
hold on
y = [0,1.4e-3];
plot(median(WTP(:))*ones(2,1),y)
hold off

Welfare = WelfareCell{3};
MarginalCost = MarginalCostCell{3};
Frictions = FrictionsCell{3};

MFXCorrelations = corr([Frictions(WTP(:) >= prctile(WTP(:),50-offset) & WTP(:) <= prctile(WTP(:),50+offset)),MarginalCost(WTP(:) >= prctile(WTP(:),50-offset) & WTP(:) <= prctile(WTP(:),50+offset)),Welfare(WTP(:) >= prctile(WTP(:),50-offset) & WTP(:) <= prctile(WTP(:),50+offset))]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Section 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define survey vars
SvyDummies = cat(5, BenefitsKnowledgeAnyIncorrect,BenefitsKnowledgeAnyNotSure,TimeCost,TimeCost.*TimeCostAccept, TimeCost.*TimeCostDontLike,ProviderNetworkHSPMore, ProviderNetworkPPOMore, ProviderNetworkNotSure,TMEGuessOverestimate, TMEGuessUnderestimate, TMEGuessNotSure,TaxBenefitsMisunderstands,TaxBenefitsNotSure);

% Desubsidize the HSP
MSsubsidy = zeros(size(UtilityHSP));
index1 = find(familysizes==1);
index2 = find(familysizes==2);
index3 = find(min(familysizes,3)==3);
MSsubsidy(index1) = 1500;
MSsubsidy(index2) = 3000;
MSsubsidy(index3) = 3750;

HSP = UtilityHSP - MSsubsidy;

% People who chose PPO
finder = find(choices(:,2)==0);
varlist = {'RA','E','HSP','UtilityPPO','ages','incomes','plans','genders','familysizes','famdraws','famdrawsOOP'};
for i = 1:length(varlist)
    eval([varlist{i},' = ',varlist{i},'(finder,2,:,:);'])
end
SvyDummies = SvyDummies(finder,2,:,:,:);

% Use parameters from Handel and Kolstad (2015) "Health Insurance for "Humans"" - Full Model with Frictions and Inertia
Params = FullModelwHSAAlpha;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%  Figures  %%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Use tier 3 for analysis
finder = find(min(familysizes(:,1,1,1),3)==3);
alphalist = 0:0.1:1;

% WTP to buy the PPO above the HSP
WTPs = cell(numel(alphalist),1);
for alphaIND = 1:numel(alphalist)
	alpha = alphalist(alphaIND);
	% Calculate WTP under adjusted frictions
	NewParams = Params;
	NewParams(7:numel(NewParams)) = (1-alpha) * NewParams(7:numel(NewParams));
	ComputedWTP = 0 - CertaintyEquivalence(NewParams,RA,E,HSP,ages,genders,incomes,familysizes,plans,SvyDummies,[],'true');
	WTPs{alphaIND} = ComputedWTP(finder,:,:);
end
OOP = famdrawsOOP(finder,:,:,:);
allowed = famdraws(finder,:,:,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Figures: Mkt Equilibrium with different friction levels %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for alpha = [0,0.5,1]
	alphaIND = find(alphalist==alpha);
	if alpha==1
		alphaSTR = '10';
	elseif alpha==0.5
		alphaSTR = '05';
	else
		alphaSTR = '00';
	end
	
	WTP = WTPs{alphaIND};
	FriclessWTP = WTPs{find(alphalist==1)};
	
	% Demand = WTP for PPO
	[Demand,ind] = sort(WTP(:),'descend');

	% Marginal cost = expected cost of the next buyer, where next is ordered by WTP.
	meanOOP = mean(OOP,3);
	MarginalCost = meanOOP(ind);
	meanCost = mean(allowed,3);
	ExpectedCost = meanCost(ind);
	
    % In HHW framework, average cost is difference in average costs between two plans for a given share
	AverageCostPPO = cumsum(ExpectedCost) ./ (1:length(ExpectedCost))';
	AverageCostHDHP = (sum(ExpectedCost - MarginalCost) - cumsum(ExpectedCost - MarginalCost)) ./ (length(ExpectedCost):-1:1)';
	AverageCost = AverageCostPPO - AverageCostHDHP;

    % Value is demand without frictions. We will plot this conditional on the employee's position on the with-frictions demand curve.
    Value = FriclessWTP(ind);
    
    %Make x-axis start at 0 and end at 1
	xaxis = (1:length(MarginalCost))' / length(MarginalCost);
	
            %%%%%%%%%%%%%%%%%%%%%%%%
            %%%%% Spline plot %%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%
	
	% Define the number of splines, then divide the x-axis accordingly
	numsplines = 10;
    
    % 11 groups, where the first and last are half the size of the other 9 groups, and form the intercept and the end point
	splinexaxis = floor((xaxis*numsplines) + .5 - 10^-6) + 1;
	
	% Mean spending/demand within each group
	DemandSpline = accumarray(splinexaxis,Demand,[],@mean);
	AverageCostSpline = accumarray(splinexaxis,AverageCost,[],@mean);
	MarginalCostSpline = accumarray(splinexaxis,MarginalCost,[],@mean);
    
    % Compute avg. value premium above MC for each group, and add it to MC
	RiskPremium = accumarray(splinexaxis,Value - MarginalCost,[],@mean);
	ValueSpline = MarginalCostSpline + RiskPremium;
	splinexaxis = ((accumarray(splinexaxis,splinexaxis,[],@max)) - 1) / numsplines;
	if alpha~=1
		plot(splinexaxis,DemandSpline,'k-',splinexaxis,AverageCostSpline,'k--',splinexaxis,MarginalCostSpline,'k-.',splinexaxis,ValueSpline,'k:')
		legend('Demand','Average Cost','Marginal Cost','Value')
	else
		plot(splinexaxis,DemandSpline,'k-',splinexaxis,AverageCostSpline,'k--',splinexaxis,MarginalCostSpline,'k-.')
		legend('Demand','Average Cost','Marginal Cost')
	end
	xlabel('Quantity')
	ylabel('Value')
	outfile = [MktEquilibriumAlpha',alphaSTR,'_HHW.png'];
	saveas(figure(1),outfile)
	eval(['Alpha',num2str(alpha*10),'AverageCostSpline = AverageCostSpline;'])
end
	
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% Figure: Avg Cost curves with different friction levels %%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot(splinexaxis,Alpha0AverageCostSpline,'k-',splinexaxis,Alpha5AverageCostSpline,'k--',splinexaxis,Alpha10AverageCostSpline,'k-.')
legend('\alpha = 0','\alpha = 0.5','\alpha = 1')

xlabel('Quantity')
ylabel('Cost')

outfile = [AvgCostCurves_HHW.png'];
saveas(figure(1),outfile)

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% Figures: Mkt equilibrium with risk adjustment %%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for alpha = [0,1]
	if alpha==0
		alphaSTR = '00';
	else
		alphaSTR = '10';
    end	
	WTP = WTPs{find(alphalist==alpha)};
	
	% Get Demand = WTP for PPO, sort descending
	[Demand,ind] = sort(WTP(:),'descend');
	AverageCostRiskAdjusted = cell(3,1);
	for betaIND = 1:3
		beta = (betaIND - 1)/2;
		
		% Marginal cost = expected cost of the next buyer, where next is ordered by WTP.
		meanOOP = (1-beta)*mean(OOP,3) + beta*mean(OOP(:));
		MarginalCost = meanOOP(ind);		
		meanCost = mean(allowed,3);
		ExpectedCost = meanCost(ind);
		
        % In HHW framework, average cost is difference in average costs between two plans for a given share
		AverageCostPPO = cumsum(ExpectedCost) ./ (1:length(ExpectedCost))';
		AverageCostHDHP = (sum(ExpectedCost - MarginalCost) - cumsum(ExpectedCost - MarginalCost)) ./ (length(ExpectedCost):-1:1)';
		AverageCostPPOPopulation = sum(ExpectedCost) / length(ExpectedCost);
		AverageCostHDHPPopulation = sum(ExpectedCost - MarginalCost) / length(ExpectedCost);
		AverageCostRiskAdjusted{betaIND} = (1-beta)*(AverageCostPPO - AverageCostHDHP) + beta*(AverageCostPPOPopulation - AverageCostHDHPPopulation);
	end
	
		splinexaxis = floor((xaxis*numsplines) + .5 - 10^-6) + 1;
	
	% Get the mean spending/demand within each group
	DemandSpline = accumarray(splinexaxis,Demand,[],@mean);
	AverageCostRiskAdjustedSpline = cell(3,1);
	for betaIND = 1:3
		AverageCostRiskAdjustedSpline{betaIND} = accumarray(splinexaxis,AverageCostRiskAdjusted{betaIND},[],@mean);
    end
    
	splinexaxis = ((accumarray(splinexaxis,splinexaxis,[],@max)) - 1) / numsplines;
	plot(splinexaxis,DemandSpline,'k-',splinexaxis,AverageCostRiskAdjustedSpline{1},'k--',splinexaxis,AverageCostRiskAdjustedSpline{2},'k-.',splinexaxis,AverageCostRiskAdjustedSpline{3},'k:')
	legend('Demand','Average Cost, \beta = 0','Average Cost, \beta = 0.5','Average Cost, \beta = 1')

	outfile = ['RiskAdjustmentMktEquilibriumAlpha',alphaSTR,'_HHW.png'];
	saveas(figure(1),outfile)
end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% Figures: 3D plots of outcomes under different policies: PPO share, premium, welfare %%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Shares = NaN(11,11);
Premiums = NaN(11,11);
Welfares = NaN(11,11);

for alphaIND = 1:11
for betaIND = 1:11
	alpha = (alphaIND - 1) / 10;	
	beta = (betaIND - 1) / 10;
	
	WTP = WTPs{alphaIND};
	
    % Get Demand = WTP for PPO, sort descending
	[Demand,ind] = sort(WTP(:),'descend');

	% Marginal cost = expected cost of the next buyer, where next is ordered by WTP.
	meanOOP = (1-beta)*mean(OOP,3) + beta*mean(OOP(:));
	MarginalCost = meanOOP(ind);
	
	AverageCostPPO = cumsum(ExpectedCost) ./ (1:length(ExpectedCost))';
	AverageCostHDHP = (sum(ExpectedCost - MarginalCost) - cumsum(ExpectedCost - MarginalCost)) ./ (length(ExpectedCost):-1:1)';
	AverageCostPPOPopulation = sum(ExpectedCost) / length(ExpectedCost);
	AverageCostHDHPPopulation = sum(ExpectedCost - MarginalCost) / length(ExpectedCost);
	AverageCost = (1-beta)*(AverageCostPPO - AverageCostHDHP) + beta*(AverageCostPPOPopulation - AverageCostHDHPPopulation);
	
	%Calculate unadjusted MC for use in welfare
	meanOOPTrue = mean(OOP,3) ;
	MarginalCostTrue = meanOOPTrue(ind);
	
	Value = WTPs{find(alphalist==1)};
	Value = Value(ind);
	
	Share = mean(Demand >= AverageCost);
	
	if Share ~=0
		Premium = min(AverageCost(Demand >= AverageCost));
	else
		Premium = Inf;
	end
	Welfare = mean((Value - MarginalCostTrue).*(Demand >= AverageCost));
	
	Shares(alphaIND,betaIND) = Share;
	Premiums(alphaIND,betaIND) = Premium;
	Welfares(alphaIND,betaIND) = Welfare;
end
end

% Calculate welfare impact as relative to (0,0)
WelfareImpacts = Welfares - Welfares(1,1);

% 3D plots
surf([0:.1:1],[0:.1:1],Shares')
colormap(gray)
xlabel('Alpha')
ylabel('Beta')
zlabel('Share')
saveas(figure(1),'3DGridShares_HHW.png')
surf([0:.1:1],[0:.1:1],Premiums')
colormap(gray)
xlabel('Alpha')
ylabel('Beta')
zlabel('Premium Difference')
saveas(figure(1),'3DGridPremiums_HHW.png')
surf([0:.1:1],[0:.1:1],WelfareImpacts')
colormap(gray)
xlabel('Alpha')
ylabel('Beta')
zlabel('Welfare Impact')
saveas(figure(1),'3DGridWelfareImpacts_HHW.png')

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%% Add. stat in different policy equilibria %%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% How many people are making mistakes?
UnderInsured = NaN(11,11);
OverInsured = NaN(11,11);
MeanCS = NaN(11,11);
StDevCS = NaN(11,11);
MeanCSBuyers = NaN(11,11);
StDevCSBuyers = NaN(11,11);
MeanCSNonbuyers = NaN(11,11);
StDevCSNonbuyers = NaN(11,11);
MeanWelfareBuyers = NaN(11,11);
StDevWelfareBuyers = NaN(11,11);
MeanWelfareNonbuyers = NaN(11,11);
StDevWelfareNonbuyers = NaN(11,11);

for alphaIND = 1:11
for betaIND = 1:11
	alphaSTR = num2str(alphaIND - 1);
	betaSTR = num2str(betaIND - 1);
	if alphaIND < 11
		alphaSTR = ['0',alphaSTR];
	end
	if betaIND < 11
		betaSTR = ['0',betaSTR];
    end
	Price = Premiums(alphaIND,betaIND);
	WTP = WTPs{alphaIND};
	WTP = WTP(:);
	Value = WTPs{find(alphalist==1)};
	Value = Value(:);
	MarginalCost = (1-beta)*mean(OOP,3) + beta*mean(OOP(:));
	MarginalCost = MarginalCost(:);
	UnderInsured(alphaIND,betaIND) = mean((WTP < Price).*(Value >= Price));
	OverInsured(alphaIND,betaIND) = mean((WTP >= Price).*(Value <= Price));
	MeanCS(alphaIND,betaIND) = mean(Value - Price);
	StDevCS(alphaIND,betaIND) = std(Value - Price);
	ksdensity(Value - Price)
	xlabel('Consumer Surplus')
	ylabel('Density')
	saveas(figure(1),['PolicyEquilibria\DensityCSAlpha',alphaSTR,'Beta',betaSTR,'_HHW.png'])
	
	if ~isinf(Price)
		MeanCSBuyers(alphaIND,betaIND) = mean(Value(WTP >= Price) - Price);
		StDevCSBuyers(alphaIND,betaIND) = std(Value(WTP >= Price) - Price);
		ksdensity(Value(WTP >= Price) - Price)
		xlabel('Consumer Surplus')
		ylabel('Density')
		saveas(figure(1),['DensityCSBuyersAlpha',alphaSTR,'Beta',betaSTR,'_HHW.png'])
	else
		MeanCSBuyers(alphaIND,betaIND) = NaN;
		StDevCSBuyers(alphaIND,betaIND) = NaN;
	end
	
	MeanCSNonbuyers(alphaIND,betaIND) = mean(Value(WTP < Price) - Price);
	StDevCSNonbuyers(alphaIND,betaIND) = std(Value(WTP < Price) - Price);
	ksdensity(Value(WTP < Price) - Price)
	xlabel('Consumer Surplus')
	ylabel('Density')
	saveas(figure(1),['DensityCSNonbuyersAlpha',alphaSTR,'Beta',betaSTR,'_HHW.png'])
	
	if ~isinf(Price)
		MeanWelfareBuyers(alphaIND,betaIND) = mean(Value(WTP >= Price) - MarginalCost(WTP >= Price));
		StDevWelfareBuyers(alphaIND,betaIND) = std(Value(WTP >= Price) - MarginalCost(WTP >= Price));
		ksdensity(Value(WTP >= Price) - MarginalCost(WTP >= Price))
		xlabel('Welfare')
		ylabel('Density')
		saveas(figure(1),['DensityWelfareBuyersAlpha',alphaSTR,'Beta',betaSTR,'_HHW.png'])

		MeanWelfareNonbuyers(alphaIND,betaIND) = mean(Value(WTP < Price) - MarginalCost(WTP < Price));
		StDevWelfareNonbuyers(alphaIND,betaIND) = std(Value(WTP < Price) - MarginalCost(WTP < Price));
		ksdensity(Value(WTP < Price) - MarginalCost(WTP < Price))
		xlabel('Welfare')
		ylabel('Density')
		saveas(figure(1),['DensityWelfareNonbuyersAlpha',alphaSTR,'Beta',betaSTR,'_HHW.png'])
	else
		MeanWelfareBuyers(alphaIND,betaIND) = NaN;
		StDevWelfareBuyers(alphaIND,betaIND) = NaN;
		MeanWelfareNonbuyers(alphaIND,betaIND) = NaN;
		StDevWelfareNonbuyers(alphaIND,betaIND) = NaN;
	end
end
end